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1. Introduction 

Recent advances in ultra-fast spectroscopy allow us to monitor the dynamics of electrons 
on a femtosecond scale. This is especially interesting for strongly correlated materials, 
such as high-temperature superconductors [H El E] , since in their case the spectroscopic 
probe is able to investigate the intra-electronic redistribution of excitation energies 
before the relaxation via the lattice starts. From the theoretical point of view, this 
is obviously a challenging problem since it requires a method capable of treating the 
relaxation dynamics of a strongly correlated system out of equilibrium. In this regard, 
a state-of-the art approach is the dynamical mean-field theory (DMFT) which has 
recently been applied [I] to the single-band Hubbard model in order to study the double- 
occupancy relaxation after laser excitation. However, for the application to real systems 
this method is rather demanding from a numerical point of view since it requires the 
self-consistent solution of complex single-impurity models driven out of equilibrium. 

In this regard, the time-dependent Gutzwiller approximation (TDGA) is a 
promising alternative since it joins the numerical simplicity of standard random phase 
approximation (RPA) with the ability to capture important many-particle effects, as 
the Mott-Hubbard transition. In a series of papers (5j EJ [8], we have developed 
the TDGA which is based on a variational Ansatz for the Hubbard model [91 [10] 
evaluated in the limit of infinite spatial dimensions pT|. This approach has recently 
been generalised for the study of multi-band Hubbard models [121 D2] and is based 
on the expansion of the Gutzwiller energy functionals which depend on the density 
matrix and variational parameters related to the atomic eigenstates. In previous work 
[HI [6] the latter have been eliminated by assuming that they instantaneously adjust to 
the density fluctuations ('antiadiabaticity approximation'). As a result, one obtains an 
energy functional which only depends on the density matrix and therefore the RPA for 
density- dependent forces [HI [15] can be applied in order to evaluate response functions. 
This (approximate) version of the TDGA has been applied successfully to the evaluation 
of dynamical correlation functions in cuprate superconductors [IE] including the optical 
conductivity [17] and the magnetic susceptibility [HI [19]. It has also been related 
to Auger spectroscopy by calculating pair excitations in one- [8] and three-band [20] 
Hubbard models. 

More recently the TDGA was extended by Schiro and Fabrizio [211 [221 [23] towards 
the inclusion of time-dependent variational parameters. Concerning the evaluation 
of response functions this approach can supersede the 'antiadiabaticity assumption', 
mentioned above, since the double-occupancy dynamics follows from a time-dependent 
variational principle. However, in Refs. [211 E2] the authors focused on the study of 
quantum quenches for systems with a homogeneous ground state. In this case, the time 
dependence is captured by the double-occupancy dynamics whereas the single-particle 
density matrix is time-independent. Recent developments consider simultaneously the 
dynamics of the double occupancy and of the density matrix [2H [25] . 

In this paper we will re-derive the TDGA for a time-dependent Gutzwiller 
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variational wave-function applied to multi-band Hubbard models. Our resulting 
equations of motion will explicitly capture the coupling between the time-dependent 
variational parameters and the density matrix. We will analyse these equations in the 
small-amplitude (i.e., linear response) limit and apply the theory to the evaluation of 
dynamical charge correlations in the single-band case. It turns out that the previous 
formulation of the TDGA j5j E] is recovered in the low-frequency limit. However, the 
incorporation of fluctuations in the time-dependent density of double occupied states 
( "doublons" ) leads to additional spectral weight above the band-like excitations in very 
good agreement with exact diagonalisation and DMFT. 

The paper is organised as follows: In Sec. 12.11 we introduce the time-dependent 
variational principle which is underlying the present work. Since the corresponding 
expectation values are evaluated with multi-band Gutzwiller wave-functions the latter 
are presented in Sec. 12.21 The evaluation of time- dependent matrix elements 
is performed in Sec. 12.31 which allows for the derivation of the Lagrangian and 
corresponding equations of motion in Sec. 12.41 Our investigations are specified for 
the single-band Hubbard model in Sec. [31 where we also discuss a two-site example 
which can be treated analytically. Finally, the small amplitude limit of the TDGA is 
derived in Sec. H]and discussed in the context of response functions in Sec.[5j Numerical 
results for the dynamical charge susceptibility are presented in Sec. [6] and compared with 
dynamical mean-field theory (DMFT) and exact diagonalisation. We finally conclude 
our investigations in Sec. [71 

2. The time-dependent Gutzwiller theory 

2.1. Variational Principle 

The time-dependent Schrodinger equation for a general time-dependent Hamiltonian 
H{t) {h= 1), can be obtained by requesting that the action 



is stationary with respect to variations of the wave function. It is in general convenient 
to perform this variation based on a real Lagrangian [26J 




which leads to equations of motion that are independent of the phase and the norm of 
the wave-functions. 

If one restricts the wave-function \^f?(zi(t))) to a certain trial form, depending on a 
set of (in general complex) 'variational parameters' Zi(t), the Euler-Lagrange equations 
for Zi(t) provide an approximation for the time evolution of the system. For example, 
restricting the wave-function to Slater determinants and using the amplitude of the 
single particle orbitals as variational parameters yields the time- dependent Hartree- 
Fock approximation. In this work, we will consider variational wave functions of the 
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Gutzwiller form [HI QUI EEE] for general multi-band Hubbard models leading to the time- 
dependent Gutzwiller approximation. 

2.2. Gutzwiller energy junctionals for multi-band Hubbard models 

We first recall some results of the conventional Gutzwiller approximation for the ground 
state properties of multi-band systems. We aim to study the physics of the following 
family of models, 

where t°j denotes the 'hopping parameters' and the operators cp a annihilate (create) 
an electron with spin-orbital index a on a lattice site i. The local Hamiltonian 

Hi-i oc = £ i;a-i, 0-2 ^o-i^i,^ (4) 

<T1,<T2 

+ C jfl,0S,eS.04gt gt g. g. 

Cl ,CT2,cr 3 ,cr 4 

is determined by the orbital-dependent on-site energies s^ auC% and by the two-particle 
Coulomb interaction ^f 71 ' " 2 ' " 3 ' " 4 . Upon introducing the eigenstates |r^ and eigenvalues 
Ei-v of (jlj) (which can be readily calculated by means of standard numerical techniques) 
the local Hamiltonian can be written as 

Hi**= X>;r|rWr'|. (5) 
r 

Multi-band Gutzwiller wave-functions have the form 

|*G)=A3|*S> = n£l*S>, (6) 

i 

where |^s) is a normalised Slater determinant and the local Gutzwiller correlator is 
defined as 

A = ^V,r'|r>«<r , l . (7) 

r,r' 

Here we introduced the variational-parameter matrix Aj ; r,r' which allows us to optimise 
the occupation and the form of the eigenstates of P;. 

The evaluation of expectations values with respect to the wave function fl6]) is 
a difficult many-particle problem, which cannot be solved in general. As shown 
in Refs. [271 EH], one can derive analytical expressions for the variational ground-state 
energy in the limit of infinite spatial dimensions (D — > oo). Using this energy functional 
for the study of finite-dimensional systems is usually denoted as the 'Gutzwiller 
approximation'. This approach is the basis of most applications of Gutzwiller wave 
functions in studies of real materials and it will also be used in this work. One 
should keep in mind, however, that the Gutzwiller approximation has its limitations 
and the study of some phenomena requires an evaluation of expectation values in finite 
dimensions [291 ED] ■ 
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For the evaluation of Gutzwiller wave functions in infinite dimensions it is most 
convenient to impose the following (local) constraints J27J 12H] 

(Ptp)* s = 1 , (8) 

(PP&l c CT ,)* s = (c+c^s • (9) 

With these constraints, the expectation value of the local Hamiltonian (J5J) reads 

(i7ioc)* G = - E ' r ^r,ri^r,r2 m ri,r2 , (10) 

r,ri,r 2 

where 

m ri ,r 2 = ((|ri>(T 2 |)>» B (11) 

can be calculated by means of Wick's theorem. 

The expectation value of a hopping operator in infinite dimensions has the form 

<4AJ*c = E «3 (&Y <€iW* s . (12) 

where the (local) renormalisation matrix (f a is an analytic function of the variational 
parameters and of the non-interacting local density matrix 

CW' = (c[a\a'h s • (13) 

The explicit form of the renormalisation matrix is given, e.g., in Ref. [31] but will not 
be needed for our further considerations in this work. In the following, we assume that 
the correlated and the non-correlated local density matrix are equal, 

C« C ;cr,cr' — (\a%a')^G = Ci- a ,a' ■ (14) 

This is the case when all non-degenerate orbitals on a lattice site belong to different 
representations of its point symmetry group. 

In summary, the expectation value of the Hamiltonian Q, 

E GA = E GA (\ ( *\p) , (15) 

is a function of all variational parameters = {A-*p r ,} and of the non-interacting 
density matrix p with the elements 

The same holds for the constraints (jHJ), ©, for which we will use the abbreviation 

g n (X { * ) ,f>)) = , l<n<n c (17) 

where n c is the maximum number of independent constraints. In Sec. [3] we apply these 
results to the special case of a single-band Hubbard model. 
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2.3. Evaluation of time- dependent matrix elements 



In this section, we will apply the concept, introduced in Sec. 12.11 to the general class of 
Gutzwiller wave functions 

i* G (t)> = p G {t)\Mt)) = n A(*)i*s(*)> , (is) 

i 

where the single-particle product states |^s(^)) an d the local correlation operators 

A(*) = X)W'(*)l r )«<r / | (19) 

r,r' 

are now time-dependent quantities. 

The state |\l/s(^)) can be written as 

l*s(f)> = W h] M ni l™> • (20) 

7 

Here, n 7 e (0, 1) determines which of the single particle states |7(i)), described by the 
operators 

^) = ^ U ^)4, (21) 

V 

are occupied and u Vjl (t) is a (time-dependent) unitary transformation, 

5Z< 1>7 (0^2,t(0 = K,v 2 ■ (22) 

7 

The functions u Vjl constitute a second set of (time-dependent) variational parameters 
(in addition to Aj ; r,r'(^)) an< l determine the single-particle wave function |^s)- Note 
that the time dependence of the operators (12T]) implies that the non-interacting density 
matrix 

Pv,v'(t) = (4'0*s(*) = ^ n 7<', 7 (*K, 7 (0 ( 23 ) 

7 

is also time dependent. 

We start with a consideration of the time derivative in Eq. (j2J) which requires the 
evaluation of 

(*g|*g) = (fr s |P G A^s) + (^s|P g ^g|^s) m 
(*g|*g) (v&s|P g Pg|^s) (*s|P g Pg|*s) 
and its complex conjugate. With Eqs. (|20|) and (121]) . we find 

|*s> = X)^Xl*s> (25) 

7 

= EE^<y^ 7 l*s). (26) 
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This equation allows us to evaluate the contribution of the first term on the r.h.s. of 
Eq. ([MD as 

(^sl^Pcl^s) 
{y s \PlP G \y s } 

^EE-x,^P^ (27) 

= ^ n^u v ^u* va . (28) 

In the last line, we have used that, in all relevant applications, Pq and |^s) have the 
same symmetry and, therefore, all contributions with 7 ^ 7' vanish in (12?]) . 

We now proceed with a consideration of the second term on the r.h.s. of 
Eq. (I24p . With the definition of the correlation operator Pq we find 

(*s|^Pg|*s> = E H ( II p i p i) A f A|*s) • (29) 

The r.h.s. of (129]) can be evaluated by the standard diagrammatic techniques in infinite 
dimensions [27J. This leads to 



- E E ^;ri,r2^;ri,r3 m i;r 2 ,r 3 , (30) 
» ri,r 2 ,r 3 

where ^i ; r 2 ,r 3 = m r,r 2 ,r 3 (t) , as defined in Eq. ( ITTj) . depends on the local elements of the 
density matrix ( 123]) . 

S.^. Lagrangian and equations of motion 

From Eqs. (1271) . (1301) . together with the expectation value of the Gutzwiller energy 
derived in Sec. I2.2[ we are now in the position to derive the Lagrangian Eq. (j2J) . However, 
we also need to include two sets of constraints, (i) the unitarity of u vr/ and (ii) the 
Gutzwiller constraints (fl7|) . Therefore we finally obtain the following Lagrangian 

L = [^*;ri,r 2 ^;ri,r 3 - A*. rijr2 Ai ;ri ,r 3 «^ ; r 2 ,r 3 (31) 

i ri,r 2 ,r 3 

+ o fcvr w *vr ~ ^.7^,7] ~~ £ GA (A W ,p) 

- J] ^7,7'W«, 7 ^Y - 1) -^A n (% n (AW,«W) 

where A n (t) and f2 7i7 /(t) are (real) Lagrange-parameters. As will be exemplified below in 
the single-band case, the original Hamiltonian can be time dependent which will reflect 
in a time dependence of E GA and allows for a coupling with arbitrary external fields. 
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From Eq. ( 1311) the Euler-Lagrange equations can now be derived in the standard 
way. The equation for the variational parameters u* v reads 

i<7 + irr- 1 ^ GA + V ^ A + yX + >J ) + >J ^,^,7 = 0(32) 



^>7 \ n J v 2 



which in terms of the density matrix f[2"3"j) can be rewritten as 

i^=[A Gfl ,p]. (33) 



ft -^( £ ° A + l/A + ^ A " 9 ") (34 » 



Here we have introduced the 'Gutzwiller Hamiltonian' 
= _? 

and a potential V" A which depends on the (time-dependent) phases of Ar,r' 

vX = 2 £ £ [^;ri,r 2 ^;ri,r 3 - ^*;ri,r 2 ^;ri,r 3 "^;r 2 ,r 3 • (35) 
i ri,r 2 ,r 3 

Note that the same equation of motion for p is obtained in the previous formulation 
of the TDGA [SI El El E]- The new ingredient in the present formulation is the 
implementation of the explicit time dependence of the variational parameters A*. r r . It 
is obtained from Eq. (I3T1) as 

i£ (^; r i- r 3 m *;r 2i r 3 + o^*; r i.r3" l i;r2,r 3 ) (36) 
r 3 ^ 



2 



3 





\ n / 



Equations (1331) and (1361) for p(t) and Aj ; r,r'(*) constitute the time-dependent Gutzwiller 
theory for multi-band Hubbard models. 

3. Example: The one-band model 

3.1. Evaluation of the time- dependent GA energy 

In order to make the results, derived in the previous section, more transparent we 
consider the case of the single band Hubbard model, 

i,j f=ti4- * 

where tij denotes the 'hopping parameters' (t^i = 0) and the operators cfl annihilate 
(create) an electron with spin index a on a lattice site i. We further introduced 

Hi,\ oc {t) = y^Vi,a(t)ni ta + Ui(t)ni^n iti , (38) 
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and hi 7a = cJ CT Ci )CT . All parameters in the Hamiltonian can be time and spatial-dependent 
which allows us to study the response to scalar fields Vi(t), vector potential fields through 
the Peierls substitution, 

= -q A-df, (39) 

and modulations of the on-site interaction. Here, we introduced q = — |e| for electrons. 

The local Hamiltonian can be diagonalised by the states \T) i = \d) i ,\a) i ,\0). 
in which the site i is double occupied (i.e. has a doublon), single occupied by an 
electron with spin a, or empty. Restricting for simplicity to paramagnetic states, we 
can work with a diagonal matrix of variational parameters, Ajrr' = Ajr^rr'- Thus the 
local 'Gutzwiller projection operator' reads, 

Pi = X d ,i\d)iM + Adt)«(tl + KH)M + A0,i|0)«(0| , (40) 

where the variational parameters A« ; r are related to the probability p^r of finding a 
configuration T at site i according to 

p tT = (^ G |r). .<r|* G ) = |A jr | 2 m ir . (41) 

We have four configuration probabilities per site which we denote as p^r = Ei, S^ a , 
Di for empty, single, and doublon occupied sites. In the present case, the constraints 
(ED,® read, 



D + S t + S x + E 



1 , (42) 



D + S a = p a 



(43) 



where 



indicates that the index i is implicit everywhere in the expression. The first 

constraint is the statement YlrPr = 1> as ^ should be, while the second constraint 
implies that local charges are the same in the correlated and the uncorrelated state. 
Obviously this also guarantees that the total charge per spin in the system is the same 
in both states. We will show below that these constrains lead to equations of motions 
that nicely respect charge conservation. 

In the real-space basis the index v in Eq. ( 1231) stands for i, a. For paramagnetic 
states the uncorrelated density matrix is diagonal with respect to the spin variables, 
Pi t a t j,a-' = Pi,j;aO~a,a' ■ We will also use the short-hand notation, p ia = Pi a ^ a . 

According to Eq. (112j) the expectation value of the one-band hopping operator in 
infinite dimensions can be written as 

(4A^)*G = '/'"'Ij.r-Pl'" ' ( 44 ) 

where the (local) renormalisation factors are given by 



q a = A* A (l - pa) + \* d \apa 



i 



(45) 
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and we used the notation 

f=| and |=| . (46) 

The q a factors renormalise the probability amplitude for the annihilation of an electron 
on site j and the creation of an electron on site i. Each one of these processes has two 
possible channels. For example the creation of an electron at i with spin up can be 
seen as a transition from an empty state to an |^) state [leading to the first term in 
Eq. (Jl5"j) ] or a transition from \l) to a doublon state [leading to second term in Eq. fj45l) ]. 
Since the variational parameters are now complex, these two channels can interfere 
either constructively or destructively affecting the total hopping amplitude. This issue 



is discussed further in |Appendix A| where the physical origin of this renormalisation is 
exemplified in a simple two-site case. 

It is convenient to write the parameters Ar in terms of a real positive amplitude 
and a phase 



A, 



— < 



,i<pr 



(47) 



which are used in the following as the dynamical variables. With these definitions the 
hopping renormalisation factors read, 



Qa 



-IX. 



-im 



with the definitions for site i, 



S a E | 
m d m^ 



Pa) 



p t -p l + D)(p a -D) 



Ar(l - Pa) 



-Pa 



Pa\ 



Pa) 



r] = ip d + y?0 - v?t - 



Xa =Wa~ <PQ 



(48) 

(49) 

(50) 

(51) 
(52) 



The phases ipi i(T , ip^i have been eliminated in favor of Xi,a and r] ija . Note that ip®^ does 
not appear anywhere in the functional and therefore can be disregarded as a dynamical 
variable. In addition, we have used the constraints Eqs. fl42T ) . fl43|) to eliminate E and 
S a in favor of D. Therefore our dynamical variables are the single-particle amplitudes 
u Vil , the double occupancy Di and the phases Xi,a and rji. 

In summary, the expectation value of the time- dependent single-band Hubbard 
model, 

(9 G (t)\H(t)\9 G (t)) 



E^ A (p, D h T] h Xi,a) 



(53) 



(tf G (f)|* G (f)> 

is a function of the variational parameters and of the non-interacting density matrix p 



E 



GA 



i,jQi,aQj,aPj,i,< 



ViaPi 



(54) 
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3.2. Lagrangian and equations of motion 

We are now in the position to evaluate the Lagrangian Eq. ( I3TT) which can be written 



as 



M°-=t,| 

u* )7 *V,7 - 1 ) (56) 



L = - ^ + Vi) - Pi,a( V i,o- + Xi,a) 

i i : a 

t i^ ix ^~ x ^\q%, h . + q d>j ,J m ){q^ + qdM^Pjw ( 55 ) 

j'0-=t,4 

w,7 v,v' \ 7 

Note that, since we have implemented the constraints (TlTl) explicitly, the corresponding 
Lagrange-parameter terms are not needed. 

The Lagrangian is invariant with respect to a gauge transformation of the form, 

XiAt) ~> XiAt) + x' i ,a( t )- 

Notice that the hopping amplitude and the site energy transform in a way that 
generalises to the lattice the usual gauge transformation in the continuum (qA — > 
qA — V%, v — >■ v + %). Hence, Xi,<r plays the role of a gauge phase and implements 
charge conservation. Indeed, the Euler-Lagrange equations for Xi,a yield the usual charge 
conservation law, 

/V = ^2kj > ( 57 ) 

3 

where the current in a bond is given by 

j id = i[t id e i{x ^~ x '' a) {q $JtU + q d j,ae irij ){q(D,i,a + qd,i,* e ~ iVi )Pi,i;v ~ hx -] ■ ( 58 ) 

The Euler-Lagrange equations for the variational parameters u* yield again the 
equation of motion for the density matrix 

ip=[h GA ,p], (59) 

with the 'Gutzwiller Hamiltonian' matrix 

> /-\,, I (00) 



and the last term ensuring gauge invariance. 

In addition to (|59|) . one obtains equations of motion for the double-occupancy 
parameters and the phases. For rji we obtain the Euler-Lagrange equation 

A = ij>,i^> «"W + Qd^^q^e^i^AjA^ ~ (61) 

JO 

From Eq. (|55p we see that rj plays for D a similar role as the gauge phase x f°r the 
charge. A time dependent r] is equivalent to a change in the Hubbard U. However, there 
are important differences. In the case of a uniform system for n > 1 and U — > oo the 
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probability to find empty sites E — > which leads to q$ t j >a ~> [c.f., Eq. (p£9]) ]. Then 
rj becomes a gauge phase and Di is conserved as can be easily checked from the charge 
constraints. Using E instead of D as variational parameter one arrives to the conclusion 
that E is conserved for n < 1 and U — > oo. In general, however, rj is not a gauge phase 
and D is of course not conserved. This reflects the fact that, when an electron jumps 
from a doubly-occupied site i to site j, one may have the process \d)i\a)j — > \a)i\d)j 
which conserves D but one can also have the process \d)i\ty)j — > \a)i\a)j which does not 
conserve D. Therefore, in general, t]i(t) is not arbitrary and has observable physical 
consequences as will be explained for the 2-site example (cf. the following section). On 
the other hand since Xi{t) is a g au g e phase we can work in a gauge in which Xi{t) = 0. 
Finally form requiring stationarity with respect to D we get 

+ h.c. (62) 
3.3. Non interacting limit 

As a check of the consistency of the equations of motion it is instructive to see how 
the non interacting limit is recovered when C/j — > 0. In this case the double occupancy 
should f actor ise as D { = p^Pii- Using this as an ansatz together with 77$ = it is easy to 
check that Eq. (!62l) is satisfied. This follows from the fact that qi$, a + qd,i,o is the hopping 
renormalisation of the static theory and attains its maximum value g©^^ + qd,i,a = 1 as 
a function of Di precisely when Di = p^p^. Thus its derivative as a function of Di 
evaluated at Di = p^Pii vanishes (as can be checked from an explicit computation) and 
Eq. (|62|) is satisfied. 

Using the same ansatz we notice that in this limit Eq. ( 16TT) can be written as, 

Di = } J qd,i,aji,j;a = P^pV (63) 
jo a 

which completes the consistency check. On the last passage, we used Eqs. f )50|) and fl57|) . 
For small U, one would recover the time dependent Hartree-Fock approximation which 
in the small amplitude limit corresponds to the usual RPA. 

3.4- Two-site example 

In order to clarify the meaning of the TDGA equations it is interesting to consider the 
following two-site, two-electron example whose exact time-dependent evolution can be 
found analytically. The Hamiltonian is assumed to be time-independent with parameters 
= 0, t\ t 2 = —to; the interaction on site 2 is infinite, ?7 2 = °°; while U% is a variable. 
Albeit simple, the problem is in the strong coupling limit and provides a non-trivial 
test for the performance of the theory. Even more, been a zero dimensional problem we 
expect it to be the more demanding test bed for the TDGA which is based on infinite 
dimension results. 
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3.4-1- Exact solution: The two-site Hamiltonian defined above is diagonalised by the 

states, 

|* ± ) = a±|d0)+6±|s) 



with 



1 

3' 



I s ) = -^( d \Ai + 4 t a u)l vac > 



and eigenvalues E± given by, 

E± = ^(U±uj ), (64) 

co = JlP + Stl - (65) 
The time-dependent wave function can be expanded as 

|tf(t)) = a + \^ + )e- iE + l + a_|^_)e- i£ - t , 

such that the double occupancy is given by 

D 1 (t) = (^(t)\S)(S\y(t)) , 

where 

(M\V(t)) = a + a + e- [E+t + a.a.e-' lE - 1 . 

Independently of the initial condition, as long as there is a finite overlap with both 
eigenstates, the double occupancy has a fluctuating part going like ~ cos(cuo^)- 
As an example we choose the starting state at t = as 

|tf(0)> = c^Jvac) = |d0) . 

The probability of finding the system in the state \d$) is then given by 

4+2 4+2 

Di(t) — 1 j[l — cos^)] ~ 1 - -f [1 - co S (cu t)] , 

u U 

where the last approximate equality is valid for U » t . Clearly the probability to find 
the system in the state \s) is 1 — D 1 from which one finds, 

m(f) = i + D 1 (t) , 

n 2 (t) = 1-D 1 (t) , 

with Hi = n it + riii- 
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Figure 1. Effective potential for charge fluctuations in the two site case. 



3.4-2. Time- dependent Gutzwiller approximation: To solve the TDGA equations it is 
useful to note that U 2 = oo leads to D 2 (t) — >• 0. If the constraints where exactly satisfied 
in the GA this would also imply E\(t) = 0. However, the numerical solution of the static 
GA shows that E\ is nonzero but very small for any U\ (Ex < 0.08) and vanishes for 
U~i — > — oo. In the following we assume for simplicity E\ — 0. Then the constraint 
Eq. (|42p implies Di — n\ — 1 and one can evaluate the TDGA equations analytically. 
The hopping renormalisation factors simplify to, 

Q2a = e _1X2CT qm : 2a , 
and the energy becomes 

E GA = -4t e- i(xiCT -* 2CT+r?l) ^1 - i J p 2 ia + h.c. + U x (m - 1) . 

The bonding single-particle state is defined by 

h\ a = u la c\ a + u 2 ac{ a . 
Without loss of generality we set u\ = x/^fe 1 ^ and u 2 = v/^- which yields 

-id, ^1(2-^1) 



and 



P21a = e 



uiu x = — + «y0 • 



The Lagrangian reads, 

L = —(Ui+rfi + 4>)ni - 4t cos(?7i + (j))v(ni) 
where we chose a gauge in which \ = and 

v(n x ) = - ( 1 - - ) y/ni(2 - n x ) . 



(66) 
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f/i/(4t) 

Figure 2. Oscillatory frequency in the TDGA for small amplitudes (dashed line) and 
in the exact solution for arbitrary fluctuations (full line) . 

Furthermore, since <fi and rji play the same role, we can set <fi = 0. Thus we get a 
Lagrangian with two dynamical variables n\ and 171 which are conjugate. Their equations 
of motion have the form, 

ri! = -4t sin(77i)t)(n 1 ) , (67) 
^ = -^-^ocos^y^x) . (68) 

The 'potential' v(ni) is shown in Fig. [TJ The static solution is given by rji = and 
v'{n\) = —Ui/(At ). Hence, for U± = the static charge is n\ = | (l + \/5) 1.62 
and it decreases towards ni = 1 when U\ — >■ t/J = 4t where a spurious Brinkman-Rice 
transition occurs. For negative Z7i instead the charge tends asymptotically to 2 when 
C/i ->■ — oo. 

Equations (1671) . (1681) can be readily solved for small oscillations around the static 
solution. We find that the oscillatory frequency is u = 4t -\A' // (n ( j ) )(— v (nf)). For 
negative t/i and until U\ = is approached we find that this frequency is in excellent 
agreement with the exact oscillation frequency of the two-site problem (c.f., Fig. [2]). 
For positive U\ as the spurious Brinkman-Rice point is approached, not surprisingly 
the approximations fail: the exact frequency increases monotonously as a function of 
U\ while the approximate one vanishes at the Brinkman-Rice point. We attribute this 
to the failure of the GA in this finite-site system. Indeed we will show below that for 
high-dimensions a similar softening approaching a Mott state is real. The softening can 
be traced back to the fact that in the Brinkman-Rice phase the doublon charge gets 
frozen at zero therefore it is conserved and 771 becomes a gauge phase which leads to an 
energy independent of r\\. We will see that a similar phenomena is found in the usual 
Brinkman-Rice transition. 
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4. Linear response: The small amplitude limit 

In the previous two-site model, the expansion of the effective potential v(ni), Eq. ( 1661) . 
around the stationary value n\ yields the dynamics in the vicinity of the GA saddle- 
point. Such an expansion is also important for the evaluation of response functions, 
where a (weak) external perturbation drives the system out of equilibrium. 

Based on our general scheme derived in Sec. I2.4[ the small amplitude limit is 
obtained by expanding the equations of motion, Eqs. (|33l) . (l36l) . around the ground 
state values p° and A° r ,, 

p(t) « p° + 5p(t) , (69) 

V,r'W~A» ir ,+dV(t). (70) 

For example, the first term on the r.h.s of Eq. (|36l) becomes 



dE GA (\(*\ p) ^ d 2 E GA C\^,~p) -, M 

— 2s 7n* n\ dX k;r 3 ,r 4 { t ) 

dX r,r u r 2 k t< A c> \r 1 ,r 2 c ' x k;r 3 ,r 4 



fc,r 3 ,r 4 

d 2 E GA (\ { *\p) 



5p v ,v'{t) ■ (71) 



After the linearisation of Eqs. ( 133|) and ( 1361) and with the harmonic Ansatz 

5X l -T,r(t) = 5X l , T , r (u)e iuJt , (72) 

SpvA*) = SPv,v'{u)e mt , (73) 

we finally end up with a linear set of equations for Ai ; r,r'(^) and 5p v>v f(u>) which can be 
solved numerically. Note that at zero frequency, the l.h.s. of Eq. ( 1361) vanishes. This 
equation then recovers exactly the 'antiadiabaticity assumption' which was used in the 
previous TDGA formulation [5j [6j [7J [8] . Within the single band Hubbard model we 
will demonstrate in the following that the inclusion of the full time-dependence of the 
variational parameters Aj ; r,r'(^) generates additional features in the dynamical charge 
correlations which are absent in the 'antiadiabatic approximation'. 

5. Response functions in the single-band Hubbard model 

In the single-band Hubbard model the three most relevant response channels are related 
to the coupling to time-dependent magnetic, charge and pair fields. This requires the 
computation of the (transversal) 'magnetic susceptibility' 

x£(t) H<4A+4At>>(*) > ( 74 ) 

the 'charge susceptibility' 

x c.(t) = ((^;n J ))(t), (75) 
and the 'pairing susceptibility' 

H<4t4;AAt>>(*) > ( 76 ) 



Linear- Response Dynamics from the Time- Dependent Gutzwiller Approximation 17 



or their respective Fourier transforms 

i f°° 

where we have introduced hi = J2 a 

5.1. The magnetic and the pairing susceptibility 

If the state \^s) is paramagnetic or is restricted to longitudinal magnetic order, the 
charge and (transverse) magnetic susceptibilities are decoupled. Moreover, if the ground 
state does not contain superconducting correlations, also the pairing fluctuations are 
decoupled from the magnetic and charge correlations. In this case, the (mixed) second 
derivative of the energy with respect to (cj.c^ ,) and vanishes. In a similar way, one 
can show that the second derivative with respect to pairing fluctuations (c\^c\ ,) and 
5\i ; d vanishes. Therefore, in both cases, the linearised differential equations (133)) and (136)) 
are decoupled and the time dependence of 5Xi-d(t) does not enter the computation of 
(transverse) magnetic and pairing correlation functions. The susceptibilities are then 
solely determined by the solution of Eq. (1331) for the single-particle density matrix and 
the present approach agrees with the previous formulation of the TDGA involving the 
'antiadiabaticity approximation' [5j El [7J [8] . Therefore we concentrate in the following on 
the investigation of the dynamical charge-charge correlation function where the present 
approach is able to capture high-energy excitations on the scale of the Hubbard repulsion 
due to the explicit time dependence of the variational parameters. 



5.2. Dynamical charge susceptibility 

As derived in Sec. [3X the time dependence of the system is governed by small deviations 
of the density matrix Sp, the double-occupancy parameters 5Di and the phase Srji from 
their saddle-point values (indicated by a '0' superscript). Note that we consider a GA 
ground state with 77° = such that Srji = r\{. The corresponding equations of motion 

9e ga 

-*=«6T- (78) 

dE GA 

SD { = — — , (79) 

OTji 

iS'p =[h GA ,5p], (80) 

have to be solved for small deviations from the GA saddle-point. For this purpose we 
expand the Gutzwiller energy functional Eq. (151)1 up to second order in the fluctuations 

E GA = E GA,0 + E GA,l + E GA,2 (gl) 

around the saddle-point value E GA '°. 
The first order contribution yields 

£ - = r r{ ^} + £(^. + ^„) (82, 
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and the derivatives have to be taken at the saddle-point. Here, the first term describes 
particle-hole excitations within the 'bare' Gutzwiller Hamiltonian. The second term 
~ 5Di vanishes due to the saddle-point condition, whereas the last term can be expressed 
in the small amplitude limit as 

dE GA Qd,i,cr Te 7 GaI V- <4 



E .0 Tlo S *-, h GA = E .0 Tie *P» , (83) 



where we have made use of Eq. flSUl) . Thus in the small amplitude limit it is convenient 
to introduce a 'displaced double occupancy' 

A = A - E n fr* (84) 

such that the dynamics of 7/j and <5A can be expressed via the second order contribution 
of E GA as 

Q E GA,2 

-rj i= °-±L^ ( 85 ) 
' d5Di 

dE GA,2 

5Di =°— . (86) 

drji 

The following evaluation of E GA ' 2 is exemplified for a homogeneous paramagnet 
but can be straightforwardly generalised to arbitrary ground states. In momentum 
space the second order expansion involves fluctuations of T] q and 5D q which are 
coupled to fluctuations of the local 5p q = J2k a ^Pk+ q ,k an d transitive 5T q = 
£*,<x [ £ l+ q + £ k] 5 Pl+ q ,k char g e densities: 

6eGA ' 2 = ^ E + v t ^ 5t : s p-« 

9 

+ E A + -4 E 

+ ^E Kq5Dq8D-. q + J £ TjVqV-q (87) 



2^ y y y 2^M 

9 9 



with 



K = W« + + 2<_ + O + + £4f>M> 



2 v JV 



«2 / r , ^9 ?2 



+ 2 * , a n L + 



°h-aS V 9 2 flg + oS J ' 



V Tp = - q \z' + z' + _) + q ° Z ' D -^, (89) 



' ^9 — 

9? = S% , (91) 



9?" = L * + K «jrt~o , (90) 

Vffl ~r (/w 
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W^E^W + W. (92) 



K q 



N 



1 

M 



- 2e°« , (93) 
L q -e q°(zl D + zl D ) + z' D (z' + z' + _)^Y^ £ k +q ( n ^) (94) 

k,cr 

where at the saddle point the renormalisation factors become site and spin independent 
(i.e. q^ i(7 = (?0, Qdia = Qd-> liu = 1°) an d the primed letters denote derivatives which 
are specified in Appendix B We have also defined the non-renormalised single-particle 
dispersion whereas the GA quasi-particle dispersion will be denoted as Sk = (gi i0 -) 2 e°. 
Then from Eqs. (1851) . (1861) the phase and double-occupancy dynamics is given by 

_ dE GA ' 2 

l N g° p b-p q + A^afK + K g 5D g , (95) 



• dE GA ' 2 1 

sb < = "Si- = m"« ' (96) 



which upon Fourier transforming in the time domain yields for the double occupancy 

1 

Here we have defined the renormalised couplings 



b, = -r=( l°'Sp, + l^iT* P°(,, W ). (97) 



^ 76? • (98) 



n T = g q (99) 
v /2Tp? 



and the 'double occupancy' Green's function 

V °^ = J%I- < 100 » 

which has poles at Q 2 = K q /M. In case of a half-filled system one obtains 

n; = 16eg[l-u a «J , (101) 

where e$ is the energy of the non-interacting system, u = U/\8eq\, and K q = 
h ^i=i cos ( a i) ■ Interestingly the Brinkman-Rice transition u — 1 can therefore be 
associated with a 'soft mode' where fl q =o — > 0. This softening is clearly associated to 
the extra gauge invariance condition that appears at the Brinkman-Rice point which 
makes the mass M to diverge. 

Equations flH7|) - fl96|) show that in the present approach the interacting electron 
problem can be mapped to an effective electron-boson problem. Electron quasiparticles 
are coupled to a boson representing fluctuations of the double occupancy ("doublon 
fluctuations") or of its conjugate variable. Notice that the mass M of the fluctuation 
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field diverges in the two situations discussed after Eq. (155]) either E — > of D — > (c.f 
Eqs. (Jl9|) . (j6T|) and fl93|) . This follows from the fact that in those cases r\ becomes a gauge 
phase and therefore gauge invariance requires that the last term of Eq. ( IHTj) vanishes. 

In analogy with electron-phonon problems the dressed fluctuations are most 
conveniently evaluated by defining a (non-interacting) susceptibility matrix 

X ee 'V) = — f dr^ ( ( l^ q(r)5p - q(0)) ( l^ q(r)<5 i^ (0)) 



.0 



iV^V4 +q + 4 (4 +q + 4) 2 ) u + e k+q -e k [ } 

and an effective interaction matrix which is composed of the bare electron-electron 
interaction and the second order 'bosonic contribution' 

[v^ o ) + (iff ) (q ' h 

The susceptibility for the interacting system is then obtained from the following RPA 
series 

X^ = ^ + ^V£x^ (103) 
Note that in the static limit oj — > the matrix Vl ee (0) is exactly the effective interaction 



obtained within the antiadiabaticity condition in Refs. [5J E] . 
6. Results 

In this section, we present results for the local dynamical charge correlation function 

^ u + E q -E -i7) 

where |0) and \q) denote ground and excited states of the single-band Hubbard model. 
We first study the low-density regime where we compare the TDGA spectra with exact 
results for the case of two particles. For higher densities we study the performance of 
the TDGA by comparing with DMFT and exact diagonalisation results. 

6.1. Low-density regime 

In the low-density limit n — > the energy of the non-interacting system is determined 
by Eq = —Bn where 2B would be the total bandwidth in case of a rectangular density 
of states. The saddle point double occupancy in this limit reads as 

n 2 1 

D° = (105) 

4(1 + *) 

which allows for the computation of the frequency of double-occupancy oscillations 
as Q q = 2B + U, i.e., it is independent of momentum. Since the coupling to these 
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fluctuations Eqs. (ESJ) , (EHJ) vanishes with n we expect the local TDGA charge correlations 
to be composed of a renormalised low-energy part for < u < 2B and a high-energy 
excitation at u = 2B + U with spectral weight proportional to the density. 

In case of two particles one can determine the eigenenergies in Eq. f |104j) from [32] 

-=-T (106) 

U E-e k -e k+q 1 ; 

where —At < e k < 4i denotes the single-particle dispersion with bandwidth 2B = 8t. The 
ground state is obtained for q = and both particles at the bottom of the band, 
i.e., E ~ —2B. A particular solution in Eq. ( I106p is obtained for Q = q = (iv,ir) at 
Eq = U so that the excited state Eq — E = 2B + U corresponds to our TDGA result 
discussed above. In addition, the exact solution involves two-particle excitations which 
are not present in the TDGA. The maximum excitation energy is obtained for q = 
which can be estimated for a rectangular density of states as u = 4B/(e l / u - 1). The 
weight of these excitations in Eq. (11061) vanishes for zero momentum transfer but clearly 
the exact solution displays high energy features in addition to the TDGA for small q 
due to the coupling between particle-hole and particle-particle excitations. 

Figure |3] displays the low-density local charge susceptibility Eq. f 1 1 6 [) . evaluated 
with TDGA and exact diagonalisation for U/t = 5 and U/t = 10, respectively. Results 
have been obtained for 2 particles on a 8 x 8 square lattice lattice (only nearest neighbor 
hopping —t). 

As anticipated above, the spectra consist of the (dominant) low-energy bandlike 
particle-hole excitations in the range < u < 2B (cf. main panels) and a high-energy 
part at u> ~ 2B + U which is resolved in the upper insets. Similar to the previous 
investigations, which were based on the 'antiadiabaticity assumption' [HE], the TDGA 
gives a very good account of the low-energy part with respect to both, energy and 
spectral weight of the excitations. The new feature, which was previously missing (5J |6] 
in the TDGA, is the high-energy part due to the explicit consideration of the double- 
occupancy time dependence. In order to estimate the associated spectral weight, we 
show in the lower insets to Fig. [3] the first moment 



Mi(w) = / duuxioM , (107) 
Jo 

which fulfills the sum rule 

Mi(oo) = — (T) . (108) 

Here (T) denotes the kinetic energy which in case of the TDGA has to be evaluated on 
the basis of the GA. From Fig.[3]it turns out that the onset of the high-energy excitations 
is accurately captured by the TDGA, however, it overestimates the associated spectral 
weight. On the other hand, this 'additional' weight is partially compensated for in 
the bandlike excitations such that the kinetic energy of the Gutzwiller approximation 
(i.e. M^oo)) again gives a good account of the exact result. 
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Figure 3. Local charge susceptibility for 2 particles on a 8 x 8 lattice with (a) U/t = 5 
and (b) U/t = 10. We compare spectra for the exact result (black, dashed) with those 
obtained within the TDGA (solid, circles). The upper insets show the high-energy 
part of the spectra. Broadening of the excitations is e = 0.02i (main panels) and 
e = 10~ 4 t (upper insets). The lower insets depict the high-energy part of the first 
moments M\(uj). 



6.2. Density dependence 

We proceed by studying the doping dependence of Xioc(^) as a function of doping which 
is shown in Fig. H]for U/t = 8 for a square lattice. The spectra again separate into low- 
energy band-like excitations and a high-energy part due to the double-occupancy time 
dependence. Starting form the low density limit the overall weight of the high-energy 




Figure 4. TDGA local charge susceptibility for a two-dimensional square lattice 
Hubbard model (U/t = 8, nearest-neighbor hopping —t) at different densities. The 
inset displays the high-energy part on an enlarged scale whereas the upper inset shows 
the imaginary part of the double-occupancy propagator. 



excitations increases approaching half filling. In addition the high-energy feature shifts 
to lower frequencies upon doping as shown in the lower inset to Fig. (j3j). 

We show also in the upper inset the imaginary part of the local double-occupancy 
propagator, i.e. 

ImLCH = Iml J>°(g,u,) (109) 

and T>(q, lo) has been defined in Eq. (11001) . 

The double-occupancy excitations evolve from Q = 2B + U (= 16t for the present 
parameters) in the limit n — > to the frequencies Q q defined in Eq. fllOip for the 
case of half-filling. In order to understand the doping dependence of the high-energy 
feature one has also to take into account the couplings Eqs. (I98l) . (l99~l) . The coupling to 
the local fluctuations, 7^ p [Eq. ( 198]) ] is continuously decreasing with the charge carrier 
concentration and dominates at all dopings except close to half-filling where it vanishes. 
On the other hand, the coupling to transitive fluctuations 7^ T [Eq. f )99|) ] is significantly 
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smaller and only weakly doping dependent. Since at half-filling the coupling between 
local and transitive fluctuations vanishes (V Tp = 0), the local charge correlations are 
only renormalised by the static density-density interactions V q . Thus at exactly half- 
filling the coupling to the double-occupancy fluctuations vanishes and the n — 1 curve 
in Fig. (B| corresponds to the 'antiadiabaticity' result derived in Refs. [5j E]. With 
decreasing doping the increasing coupling between local density and double fluctuations 
increases and induces the shift of the high-energy feature to large frequencies. 

In order to check the quality of the TDGA at larger doping vs. other approaches we 
compare our results with dynamical mean-field theory (DMFT) [33]. Despite freezing the 
spatial fluctuations beyond mean-field, DMFT takes fully into account the local quantum 
dynamics and it is in particular reliable to describe the evolution of the spectral weight 
as a function of temperature and other control parameters like doping. DMFT maps 
the lattice model onto an Anderson Impurity model, which is solved with an 'impurity 
solver', which in the present work is exact diagonalisation [31]. Within this method the 
bath of the AIM is discretised into Nb levels, which here is taken to be 9. A Matsubara 
grid defined by an effective temperature (3 = 80 /t is used and the stability of the results 
as a function of the two control parameters has been checked. Within DMFT, the local 
dynamical correlation functions can be obtained without further approximation. As a 
result of the discrete bath, the spectra appear more spiky than in the actual solution, but 
it has been shown that this approach describes accurately the evolution of the spectral 
weight (for instance of the optical conductivity) as a function of the various control 
parameters [351 EE] . 



n=0.6 




co/t 



Figure 5. DMFT local charge susceptibility for a two-dimensional square lattice 
Hubbard model (U/t — 8, nearest-neighbor hopping —t) at densities n — 0.6, n — 0.8, 
and n = 1. 
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Figure El (main panel) shows the local charge susceptibility obtained within DMFT 
for different fillings and U/t = 8. Despite the peaky structure (which hampers a direct 
comparison with the TDGA results) it is obvious that the main features correspond 
to those obtained within the TDGA, i.e., bandlike excitations on the scale of 8t 
and additional higher energy excitations which soften and gain spectral weight upon 
increasing density and approaching the insulating phase. 

In order to test the performance of the TDGA we show in Fig. El (panels b,d) a 
comparison of the first moment Mi (a;), Eq. (I107p . evaluated within TDGA (black solid 
lines) and DMFT (green circles) for U/t = 8 and densities n = 0.6 (panel b) and n = 1.0 
(panel d). As anticipated above the TDGA gives a rather good account of the spectral 
weight evolution at lower densities where it is in good agreement with the DMFT data 
(Fig. El panel b) given the uncertainties due to the finite number of bath states. However, 
due to the vanishing coupling between electrons and double-occupancy fluctuations at 
half-filling, all the TDGA spectral weight is contained in the band-like excitations in 
this limit so that the 'antiadiabaticity' result of Refs. [5j E] is recovered. Therefore the 
corresponding first moment increases much faster than DMFT but nevertheless both 
moments approach for u — > oo due to the agreement in the kinetic energies. One should 
note that at half-filling the GA ground state actually corresponds to a spin-density wave 
and also for such symmetry-broken states we find that at half-filling the antiadiabaticity 
result of Refs. [SI E] is valid. In fact, the TDGA charge excitations on top of such a 
ground state are in reasonable agreement with exact diagonalisation as demonstrated 
in Fig. 5 of Ref . [6] . 

For comparison, we also show in Fig. [6] the result of HF+RPA theory together with 
the spectra of our previous 'approximate' TDGA [SI E] where the double-occupancy 
fluctuations have been antiadiabatically eliminated. As discussed above, at half-filling 
the antiadiabatic approximation agrees with the 'exact' evaluation of the TDGA (cf. 
panels c,d in Fig. E])- The difference becomes pronounced at lower doping where the high- 
energy feature gets significant spectral weight and due to repulsion induces a softening of 
the bandlike excitations. From panel a of Fig. El it is clear that both effects are essential 
for reproducing the very good agreement of the first moment with the DMFT result at 
n = 0.6 whereas the approximate TDGA interpolates the spectral weight between the 
bandlike and high-energy excitations. Note, however, that for small frequencies u> — > 
the approximate result agrees with the 'exact' TDGA as has already been discussed in 
Sec. HI In addition, the first moment of both approaches agrees in the limit u — > oo 
since it is set by the static GA kinetic energy. 

In case of the conventional HF+RPA approach, where the local charge susceptibility 
is given by 

HF+RPA/ ,\ _ 1 n 1 „x 

and Xqi u ) nas the same structure than the (ll)-element of Eq. (I102I) but with non- 
renormalised single-particle dispersions e° k . Since HF+RPA obviously does not capture 
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Figure 6. Imaginary part of the local charge susceptibility at densities n = 0.6 
(panel a) and n — 1 (panel c) evaluated with the TDGA (solid black), TDGA with 
antiadiabaticity (dashed red), and HF+RPA (dotted blue). Panels c,d report the 
corresponding evolution of the first moment Eq. (|107[) where also the results from 
DMFT (circles, green) are shown. Onsite repulsion: U/t = 8. 



the double-occupancy dynamics, the local charge fluctuations originate from the band- 
type excitations which are renormalised due to the RPA denominator and high frequency 
excitations are absent. Moreover, HF overestimates the kinetic energy so that the first 
moment of Xio! +RPA overshoots both the TDGA and DMFT result. As can be seen from 
Fig. M this discrepancy is most apparent close to half-filling where the renormalisation 
of the kinetic energy due to correlation effects is more pronounced. Upon reducing the 
band filling the low-energy part of the HF+RPA spectra approaches the TDGA result 
where, however, the latter approach additionally captures the high frequency excitations 
at 2B + U with spectral weight ~ n. 



7. Conclusions 



In this paper we have developed the time-dependent Gutzwiller approximation for multi- 
band Hubbard models. This approach is based on a time-dependent variational principle 
where expectation values are evaluated with the Gutzwiller variational wave-function in 
the limit of infinite dimensions. In contrast to the standard Gutzwiller approximation 
[HI HOI EH] both, variational parameters and the underlying Slater determinant, acquire 
a time dependence. In this regard our calculations generalise earlier investigations 
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by Schiro and Fabrizio [2H [22] who have studied quantum quenches in homogeneous 
systems where the time dependence of the density matrix does not couple to that of 
the variational parameters. On the other hand, momentum (or space) dependent out-of 
equilibrium displacements of the system require such a coupling as evident from our 
generalised equations of motion Eqs. (I33l).fl36lh 

We have applied this theory in the small-amplitude, i.e., linear-response limit 
and exemplified for the case of dynamical charge correlations in the single-band 
Hubbard model. In an earlier formulation of the TDGA the so-called 'antiadiabaticity 
approximation' [5j [6] has been applied, where the dynamics of the double-occupancy 
parameters was slaved by that of the density matrix. In contrast, the present approach 
explicitly incorporates the time dependence of the double-occupancy variational 
parameters which agrees with the previous formulation in the static limit. In addition it 
improves the theory in Refs. [6] by incorporating the high-energy features which are 
on the scale of the Hubbard repulsion for small densities and which position is in good 
agreement with that of exact diagonalisation. On the other hand, the spectral weight of 
the high-energy excitations is overestimated within the TDGA although it significantly 
improves the standard HF+RPA approach in this regard. Further refinement of the 
theory could be achieved by including the coupling between particle-hole and particle- 
particle excitations which have been studied in Ref. [8] in the framework of the GA. 

It is interesting that in the present approach the Brinkman-Rice transition appears 
signaled by a collective mode whose frequency goes to zero. This is not due to the 
doublon fluctuation stiffness becoming soft but because the mass of the fluctuations 
diverges. We have shown that this divergence appears each time the double occupancy 
becomes a conserved quantity which is the case in the Brinkman-Rice case where D = 0. 
It remains to be seen which of these feature remain in an exact description although 
the similarity of the TDGA results with dynamical mean field theory (DMFT) suggest 
that at least in an approximate way this physics survives in real Mott transitions. 

Within the DMFT it is quite diffult to study systems in which the momentum 
dependence of collective excitations is important as, for example, spin waves in insulators 
|37j . In such cases, the TDGA provides us with an important additional tool which 
complements the DMFT. 
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Appendix A. Physical meaning of the phases 

To understand the physical meaning of the phases appearing in Eq. (1481) and how they 
affect the hopping amplitude consider the following two-site example in which the non 
interacting state is uniform, 

l*s> = ^(4, t + 4, t )( a u + 4,|)l v ac) 

but the projectors are given by Eq. (14*01) with all Ar = 1 except \2d = —1 which 
corresponds to r] 2 = 7r and the other phases zero leading to q 2 , a = 0, i.e., destructive 
interference. The projected wave function reads, 

l*G> = \( d lA,l _ d lA,l + d lA,i + 4,t 6 U)l vac ) • 

The exact off-diagonal density matrix in the Gutzwiller wave function is given by the 
overlap between the states, 

ci,tl*o) = ^(c f i,t + 4,t)l vac ) ' 

C2,tl*e) = ^(4,t -4,t)l vac ) • 

We see that in this zero dimensional example the 'background' electron remains with the 
'wrong' phase (or momentum) leading to zero overlap, in accord with the GA derived 
in infinite dimensions. Notice, however, that if also A^ = —1 the overlap is finite in 
the exact evaluation while it is zero in the GA. Clearly such kind of process depends 
on the correlation between the phases on different sites and can not be captured by the 
factorised form ~ q* a q2,a of the GA. 

Appendix B. Derivatives of the renormalisation factors 

In Sec. 15.21 we have introduced derivatives of the hopping renormalisation factor q ij(T , 
evaluated at the saddle-point of a paramagnetic system. These are defined as follows: 

dgjg _ , dq ia _ , dq^_ _ , 
dpua ' dpu-o + ~' dDi D 
d 2 g i(7 _ d 2 q ia _ » d 2 q ia 

dpla ~ ++ ' dpuadpu-a ~ + "' dp\_ a 

d 2 q i(7 _ ,< d 2 q ia __ „ d 2 q ia 
3D 2 ~ Z °> dp^dDi ~ Z+D ' dpu^dDi 
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